
# Analysis appendix

rm(list = ls())

# This script contains the code to generate the figures and tables of the main manuscript.

#Set your path here

setwd("")

#Install and load required packages
library("pacman")
p_load(patchwork, texreg, estimatr, marginaleffects, grid, gridExtra, fixest, sf, haven, 
       ggthemes, ggdist, tidyquant, lubridate, tidyverse, corrplot, ggspatial) 

# A. Data information -------------------------------

# Figure A1

loadRData <- function(fileName){
  load(fileName)
  get(ls()[ls() != "fileName"])
}

df <- loadRData("replication/database/short_time.RData")

vars <- c("depo_ratio_te", "agriculture_workers", "services_workers", "unemployed", "var_older_te", "var_younger_te", "population_log")
data <- df[, vars]
data <- na.omit(data)
rm(df)

cor.mtest <- function(mat, ...) {
  mat <- as.matrix(mat)
  n <- ncol(mat)
  p.mat<- matrix(NA, n, n)
  diag(p.mat) <- 0
  for (i in 1:(n - 1)) {
    for (j in (i + 1):n) {
      tmp <- cor.test(mat[, i], mat[, j], ...)
      p.mat[i, j] <- p.mat[j, i] <- tmp$p.value
    }
  }
  colnames(p.mat) <- rownames(p.mat) <- colnames(mat)
  p.mat
}

data <- data |> 
  dplyr::rename("Depopulation (ratio)" = depo_ratio_te, "% workers agriculture" = agriculture_workers, "% workers services" = services_workers, "% unemployed" = unemployed, "\u0394 % 60 y/o or older" = var_older_te, "\u0394 % 16 y/o or younger" = var_younger_te, "(Log) Population" = population_log)

M <- cor(data)
p.mat <- cor.mtest(data)

png(file="replication/outcomes/images/FigureA1.png",width=1800, height=1800, pointsize=30)
col <- colorRampPalette(c("#BB4444", "#EE9988", "#FFFFFF", "#77AADD", "#4477AA"))
corrplot(M, method="color", col=col(200),  
         type="upper", order="hclust", 
         addCoef.col = "black", # Add coefficient of correlation
         tl.col="black", tl.srt=45, #Text label color and rotation
         # Combine with significance
         # hide correlation coefficient on the principal diagonal
         diag=FALSE 
)
dev.off()

# Figure A2
rm(list = ls())


loadRData <- function(fileName){
  load(fileName)
  get(ls()[ls() != "fileName"])
}

df <- loadRData("replication/database/long_time.RData")

df <- df |> 
  filter(year %in% c(1989, 1993, 1996, 2000, 2004, 2008, 2011, 2015, 2016, 2019))

sd <- df |> 
  group_by(year) %>% 
  summarise(min_sd = -(0.5*sd(depo_ratio_te, na.rm = T)),
            max_sd = 0.5*sd(depo_ratio_te, na.rm = T))

ggplot() +
  geom_histogram(data = df, aes(x = depo_ratio_te, y = after_stat(density)), bins = 100, na.rm = T, fill = "#636363") +
  scale_x_continuous(limits = c(-70, 70)) +
  geom_vline(data = sd, aes(xintercept = min_sd), color = "black") +
  geom_vline(data = sd, aes(xintercept = max_sd), color = "black") +
  facet_wrap(~year, ncol = 3) +
  theme_minimal() +
  labs(x = "Depopulation", y = "Density") +
  geom_rect(data = sd, aes(xmin = -Inf, xmax = min_sd, fill = "Increase"), colour = NA, ymin = 0, ymax = 0.045, alpha = .4) +
  geom_rect(data = sd, aes(xmin = min_sd, xmax = max_sd, fill = "No change"), colour = NA, ymin = 0, ymax = 0.045, alpha = .4) +
  geom_rect(data = sd, aes(xmin = max_sd, xmax = Inf, fill = "Decrease"), colour = NA, ymin = 0, ymax = 0.045, alpha = .4) +
  scale_fill_manual(" ",
                    values = c("#fc9272","#a1d99b", "#ffeda0"),  
                    guide = guide_legend(override.aes = list(alpha = 1))) +
  geom_text(data = sd, aes(x = min_sd - 20, y = 0.04), label = "- 0.5 SD") +
  geom_text(data = sd, aes(x = max_sd + 20, y = 0.04), label = "+ 0.5 SD") +
  theme(axis.title.x = element_text(size = 14),
        axis.title.y = element_text(size = 14),
        legend.text = element_text(size = 13),
        legend.position = "bottom", 
        strip.text = element_text(face = "bold", size = 12))

ggsave(filename = "replication/outcomes/images/FigureA2.png", height = 12, width = 10)

# Figure A3
rm(list = ls())

loadRData <- function(fileName){
  load(fileName)
  get(ls()[ls() != "fileName"])
}

df <- loadRData("replication/database/short_time.RData")

df <- df |> 
  filter(year %in% c(2011, 2015, 2016, 2019))

colour <- c("#de2d26", "#a1d99b","#3182bd")
names(colour) <- c("Cities", "Towns and Suburbs", "Rural areas")

graph <-boxplot(df$depo_ratio_te, col="skyblue", frame.plot=F) # Remove outliers
df_without <- df[!(df$depo_ratio_te %in% graph$out), ]

df_without |> 
  filter(!is.na(degurba)) |> 
  ggplot(aes(x = degurba, y = depo_ratio_te)) + 
  facet_wrap(~year) +
  ggdist::stat_halfeye(aes(fill = degurba),
                       adjust = .5, 
                       width = .6, 
                       .width = 0, 
                       justification = -.3, 
                       point_colour = NA) + 
  geom_boxplot(
    width = .25, 
    outlier.shape = 3
  ) +
  geom_point(aes(color = degurba),
             fill = "black",
             size = 1.5,
             alpha = .09,
             position = position_jitter(
               seed = 1, width = .1
             )) + 
  scale_fill_manual(values = colour) +
  scale_color_manual(values = colour) +
  geom_hline(yintercept = 0, alpha = .5, linetype = "dashed", color = "black") +
  labs(y = "% depopulation", x = " ") +
  coord_flip() +
  theme_minimal() +
  theme(axis.text.y = element_text(size = 14), 
        axis.text.x = element_text(size = 14), 
        legend.position = "none",
        strip.text = element_text(face = "bold", size = 12))

ggsave(file = "replication/outcomes/images/FigureA3.png", width = 12, height = 10)

# B. Political parties and context in Spain -------------------------------

# Figure B1
rm(list = ls())

mp <- read.csv("Data/source/cmp/MPDataset_MPDS2022a.csv")
mp <- mp[mp$countryname == "Spain", ]

mp <- mp[mp$partyname %in% c("Citizens - Party of the Citizens", "Communist Party of Spain", "People's Party", "Popular Alliance", "Popular Unity", "United Left","Spanish Socialist Workers’ Party", "United We Can", "Voice", "We can"), ]

mp <- mp[, c(5,8,168)]
mp <- mp |>  pivot_wider(values_from = rile, names_from = partyname)

mp[is.na(mp)] <- 0
mp$PP <- mp$`People's Party` + mp$`Popular Alliance`
mp$PSOE <- mp$`Spanish Socialist Workers’ Party`
mp$Vox <- mp$Voice
mp$Cs <- mp$`Citizens - Party of the Citizens`
mp$`United Left`[mp$edate %in% c("10/11/2019", "28/04/2019") ] <- 0
mp$Podemos <- mp$`Communist Party of Spain` + mp$`We can` + mp$`United Left`

mp <- mp[, c("edate","PSOE", "PP", "Vox", "Cs", "Podemos")]

mp <- mp %>% 
  pivot_longer(cols = 2:6, names_to = "party", values_to = "rile") %>% 
  separate(edate, c("day", "month", "year"), sep = "/") %>% 
  select(year, party, rile)
mp$rile[mp$rile == 0] <- NA

colour <- c("#d7301f", "#2b8cbe", "#41ae76", "#8c6bb1", "#fe9929")
names(colour) <- c("PSOE", "PP", "Vox", "Podemos", "Cs")

mp %>% 
  group_by(year, party) %>% 
  summarise(mean = mean(rile)) %>% 
  ggplot(aes(x = year, y = mean, group = party)) +
  geom_line(aes(color = party), size = 1) + 
  geom_point(aes(color = party), size = 2) +
  scale_color_manual(values = colour, name = "Party") +
  labs(x = "Year", y = "Right - left position") +
  theme(axis.title.x = element_text(size = 14),
        axis.title.y = element_text(size = 14),
        legend.text = element_text(size = 14)) +
  theme_minimal() 

ggsave(filename = "replication/outcomes/images/FigureB1.png", width = 15, height = 10)

# Figure B2
rm(list = ls())

mp <- read.csv("Data/source/cmp/MPDataset_MPDS2022a.csv")
mp <- mp[mp$countryname == "Spain", ]

mp <- mp[mp$partyname %in% c("Citizens - Party of the Citizens", "Communist Party of Spain", "People's Party", "Popular Alliance", "Popular Unity", "United Left","Spanish Socialist Workers’ Party", "United We Can", "Voice", "We can"), ]

mp <- mp[, c(5,8,171)]

mp <- mp %>% pivot_wider(values_from = welfare, names_from = partyname)

mp[is.na(mp)] <- 0
mp$PP <- mp$`People's Party` + mp$`Popular Alliance`
mp$PSOE <- mp$`Spanish Socialist Workers’ Party`
mp$Vox <- mp$Voice
mp$Cs <- mp$`Citizens - Party of the Citizens`
mp$`United Left`[mp$edate %in% c("10/11/2019", "28/04/2019") ] <- 0
mp$Podemos <- mp$`Communist Party of Spain` + mp$`We can` + mp$`United Left`

mp <- mp[, c("edate","PSOE", "PP", "Vox", "Cs", "Podemos")]

mp <- mp %>% 
  pivot_longer(cols = 2:6, names_to = "party", values_to = "welfare") %>% 
  separate(edate, c("day", "month", "year"), sep = "/") %>% 
  select(year, party, welfare)

mp$welfare[mp$welfare == 0] <- NA

colour <- c("#d7301f", "#2b8cbe", "#41ae76", "#8c6bb1", "#fe9929")
names(colour) <- c("PSOE", "PP", "Vox", "Podemos", "Cs")

mp %>% 
  group_by(year, party) %>% 
  summarise(mean = mean(welfare)) %>% 
  ggplot(aes(x = year, y = mean, group = party)) +
  geom_line(aes(color = party), linewidth = 1) + 
  geom_point(aes(color = party), size = 2) +
  scale_color_manual(values = colour, name = "Party") +
  theme_bw() +
  labs(x = "Year", y = "Welfare") +
  theme(axis.title.x = element_text(size = 14),
        axis.title.y = element_text(size = 14),
        legend.text = element_text(size = 14)) +
  theme_minimal() +
  geom_hline(yintercept = 6.5)

ggsave(filename = "replication/outcomes/images/FigureB2.png", width = 15, height = 10)

# Figure B3
rm(list = ls())

mp <- read.csv("Data/source/cmp/MPDataset_MPDS2022a.csv")
mp <- mp[mp$countryname == "Spain", ]

mp <- mp[mp$partyname %in% c("Citizens - Party of the Citizens", "Communist Party of Spain", "People's Party", "Popular Alliance", "Popular Unity", "United Left","Spanish Socialist Workers’ Party", "United We Can", "Voice", "We can"), ]

mp <- mp[, c(5,8,166)]
mp <- mp %>% pivot_wider(values_from = per703_1, names_from = partyname)

mp[is.na(mp)] <- 0
mp$PP <- mp$`People's Party` + mp$`Popular Alliance`
mp$PSOE <- mp$`Spanish Socialist Workers’ Party`
mp$Vox <- mp$Voice
mp$Cs <- mp$`Citizens - Party of the Citizens`
mp$`United Left`[mp$edate %in% c("10/11/2019", "28/04/2019") ] <- 0
mp$Podemos <- mp$`Communist Party of Spain` + mp$`We can` + mp$`United Left`

mp <- mp[, c("edate","PSOE", "PP", "Vox", "Cs", "Podemos")]

mp <- mp %>% 
  pivot_longer(cols = 2:6, names_to = "party", values_to = "agrarian") %>% 
  separate(edate, c("day", "month", "year"), sep = "/") %>% 
  select(year, party, agrarian)

mp$agrarian[mp$agrarian == 0] <- NA
mp <- na.omit(mp)

colour <- c("#d7301f", "#2b8cbe", "#41ae76", "#8c6bb1", "#fe9929")
names(colour) <- c("PSOE", "PP", "Vox", "Podemos", "Cs")

mp %>% 
  group_by(year, party) %>% 
  summarise(mean = mean(agrarian)) %>% 
  ggplot(aes(x = year, y = mean, group = party)) +
  geom_line(aes(color = party), linewidth = 1) + 
  geom_point(aes(color = party), size = 2) +
  scale_color_manual(values = colour, name = "Party") +
  theme_minimal() +
  labs(x = "Year", y = "Agrarian and farmers")

ggsave(filename = "replication/outcomes/images/FigureB3.png", width = 15, height = 10)

# Figure B4
rm(list = ls())

mp <- read.csv("Data/source/cmp/MPDataset_MPDS2022a.csv")
mp <- mp[mp$countryname == "Spain", ]

mp <- mp[mp$partyname %in% c("Citizens - Party of the Citizens", "Communist Party of Spain", "People's Party", "Popular Alliance", "Popular Unity", "United Left","Spanish Socialist Workers’ Party", "United We Can", "Voice", "We can"), ]

mp <- mp[, c(5,8,71)]

mp <- mp %>% pivot_wider(values_from = per604, names_from = partyname)

mp[is.na(mp)] <- 0

mp$PP <- mp$`People's Party` + mp$`Popular Alliance`
mp$PSOE <- mp$`Spanish Socialist Workers’ Party`
mp$Vox <- mp$Voice
mp$Cs <- mp$`Citizens - Party of the Citizens`
mp$`United Left`[mp$edate %in% c("10/11/2019", "28/04/2019") ] <- 0
mp$Podemos <- mp$`Communist Party of Spain` + mp$`We can` + mp$`United Left`

mp <- mp[, c("edate","PSOE", "PP", "Vox", "Cs", "Podemos")]

mp <- mp %>% 
  pivot_longer(cols = 2:6, names_to = "party", values_to = "traditional") %>% 
  separate(edate, c("day", "month", "year"), sep = "/") %>% 
  select(year, party, traditional)

# a$traditional[a$traditional == 0] <- NA
mp$traditional[mp$party == "Vox" & mp$year < 2019 ] <- NA
mp$traditional[mp$party == "Cs" & mp$year < 2015 ] <- NA
mp <- na.omit(mp)

colour <- c("#d7301f", "#2b8cbe", "#41ae76", "#8c6bb1", "#fe9929")
names(colour) <- c("PSOE", "PP", "Vox", "Podemos", "Cs")

mp %>% 
  group_by(year, party) %>% 
  summarise(mean = mean(traditional)) %>% 
  ggplot(aes(x = year, y = mean, group = party)) +
  geom_line(aes(color = party), size = 1) + 
  geom_point(aes(color = party), size = 2) +
  scale_color_manual(values = colour, name = "Party") +
  theme_minimal() +
  labs(x = "Year", y = "Traditional values: Negative")

ggsave(filename = "replication/outcomes/images/FigureB4.png", width = 15, height = 10)

# Figure B5
rm(list = ls())

loadRData <- function(fileName){
  #loads an RData file, and returns it
  load(fileName)
  get(ls()[ls() != "fileName"])
}

df <- loadRData("replication/database/long_time.RData")

df <- df[, c(1,3,5:11, 15)]

df <- df |> pivot_longer(cols = 3:9, names_to = "Party", values_to = "Vote")

df$Party[df$Party == "psoe" ] <- "PSOE"
df$Party[df$Party == "pp" ] <- "PP"
df$Party[df$Party == "vox" ] <- "Vox"
df$Party[df$Party == "podemos" ] <- "Podemos"
df$Party[df$Party == "cs" ] <- "Cs"
df$Party[df$Party == "panes" ] <- "NSWP"
df$Party[df$Party == "empty" ] <- "Empty"

df <- df |> filter(Party %in% c("PSOE", "PP", "Vox", "Podemos", "Cs"))

df$year <- as.character(df$year)
df$year[df$year == "2020"] <- "2019 November"
df$year[df$year == "2019"] <- "2019 April"

df$log_re <- cut(df$population_log, seq(from = 2, to = 15, by = 1))
df$log_re [df$population_log < 2 ] <- "(2,3]"

df <- df |> 
  group_by(year, Party, log_re) |> 
  summarise(mean = mean(Vote)) |> 
  drop_na()

colour <- c("#d7301f", "#2b8cbe", "#41ae76", "#8c6bb1", "#fe9929")
names(colour) <- c("PSOE", "PP", "Vox", "Podemos", "Cs")

df |> 
  ggplot(aes(x = log_re, y = mean, group = Party) ) +
  geom_line(aes(colour = Party), linewidth = 1, show.legend = T) +
  labs(x = "(Log) Population", y = "% vote") +
  facet_wrap(~ year, ncol = 3) +
  scale_color_manual(values = colour ) +
  scale_y_continuous(minor_breaks = NULL) +
  theme_minimal() +
  theme(strip.text = element_text(face = "bold", size = 18),
        legend.position = "bottom",
        legend.title = element_text(size = 16),
        legend.text = element_text(size = 16),
        axis.text.x = element_text(angle = 45, vjust = 0),
        axis.title.x = element_text(vjust = -5),
        axis.title.y = element_text(size = 16)) +
  guides(linetype = guide_legend(title= "Parties") )

ggsave(file = "replication/outcomes/images/FigureB5.png", width = 15, height = 10)

# Figure B6
rm(list = ls())

cis19n <- read_sav("replication/database/CIS/nov19.sav")
cis19a <- read_sav("replication/database/CIS/apr19.sav")
cis16 <- read_sav("replication/database/CIS/jun16.sav")
cis15 <- read_sav("replication/database/CIS/dec15.sav")
cis11 <- read_sav("replication/database/CIS/nov11.sav")

  # November 2019

    # Age
cis19n$age <- cis19n$C10
cis19n$age_re <- cut(cis19n$age, seq(from = 18, to = 83, by = 5))
cis19n$age_re[cis19n$age > 78 ] <- "(78,83]"
table(cis19n$age_re)

    # Vote
cis19n$vote <- cis19n$B22
cis19n$vote <- car::recode(cis19n$vote, "c(4)=3 ; c(5,6,7,10,20,21,22)=4; c(18) = 5")
cis19n$vote <- factor(cis19n$vote,
                      levels = c(1,2,3,4,5),
                      labels = c("PP", "PSOE", "Cs", "Podemos", "Vox"))
cis19n <- cis19n[cis19n$vote %in% c("PP", "PSOE", "Cs", "Podemos", "Vox"), ]
cis19n$year <- "Nov 2019"
cis19n <- cis19n[, c("year", "vote", "age_re")]

  # April 2019

    # Age
cis19a$age <- cis19a$P42
cis19a$age_re <- cut(cis19a$age, seq(from = 18, to = 83, by = 5))
cis19a$age_re[cis19a$age > 78 ] <- "(78,83]"
table(cis19a$age_re)

    # Vote
cis19a$vote <- cis19a$P23
cis19a$vote <- car::recode(cis19a$vote, "c(4)=3 ; c(6,7,21)=4; c(18) = 5")
cis19a$vote <- factor(cis19a$vote,
                      levels = c(1,2,3,4,5),
                      labels = c("PP", "PSOE", "Cs", "Podemos", "Vox"))
cis19a <- cis19a[cis19a$vote %in% c("PP", "PSOE", "Cs", "Podemos", "Vox"), ]
cis19a$year <- "Apr 2019"
cis19a <- cis19a[, c("year", "vote", "age_re")]

  # June 2016

    # Age
cis16$age <- cis16$P48
cis16$age_re <- cut(cis16$age, seq(from = 18, to = 83, by = 5))
cis16$age_re[cis16$age > 78 ] <- "(78,83]"
table(cis16$age_re)

# Voto
cis16$vote <- cis16$RECUERDO16
cis16$vote <- car::recode(cis16$vote, "c(3,5,6,9) = 3")
cis16$vote <- factor(cis16$vote,
                     levels = c(1,2,3,4),
                     labels = c("PP", "PSOE", "Podemos", "Cs"))
cis16 <- cis16[cis16$vote %in% c("PP", "PSOE", "Cs", "Podemos"), ]
cis16$year <- "2016"
cis16 <- cis16[, c("year", "vote", "age_re")]

  # December 2015

    # Age
cis15$age <- cis15$EDAD
cis15$age_re <- cut(cis15$age, seq(from = 18, to = 83, by = 5))
cis15$age_re[cis15$age > 78 ] <- "(78,83]"
table(cis15$age_re)

    # Vote
cis15$vote <- cis15$P31R
cis15$vote <- car::recode(cis15$vote, "c(3,5,6,7,10) = 3")
cis15$vote <- factor(cis15$vote,
                     levels = c(1,2,3,4),
                     labels = c("PP", "PSOE", "Podemos", "Cs"))
cis15 <- cis15[cis15$vote %in% c("PP", "PSOE", "Cs", "Podemos"), ]
cis15$year <- "2015"
cis15 <- cis15[, c("year", "vote", "age_re")]

  # 2011

    # Age
cis11$age <- cis11$EDAD
cis11$age_re <- cut(cis11$age, seq(from = 18, to = 83, by = 5))
cis11$age_re[cis11$age > 78 ] <- "(78,83]"
table(cis11$age_re)

  # Voto
cis11$vote <- cis11$p36
cis11$vote <- car::recode(cis11$vote, "c(3,11) = 3")
cis11$vote <- factor(cis11$vote,
                     levels = c(1,2,3),
                     labels = c("PP", "PSOE", "Podemos"))
cis11 <- cis11[cis11$vote %in% c("PP", "PSOE", "Podemos"), ]
cis11$year <- "2011"
cis11 <- cis11[, c("year", "vote", "age_re")]

  # Total
cis <- rbind(cis19n, cis19a, cis16, cis15, cis11)
rm(cis19n, cis19a, cis16, cis15, cis11)

cis <- cis |> 
  group_by(year, age_re, vote) |> 
  summarise(n = n()) |> 
  mutate(freq = (n / sum(n))*100 ) |> 
  drop_na()
cis$freq <- round(cis$freq, 2)
cis <- cis[, -4] 

mean <- cis |> 
  group_by(vote, age_re) |> 
  summarise(mean = mean(freq))

mean$year <- "Mean"
mean <- mean[, c(4,2,1,3)]
colnames(mean)[4] <- "freq"
mean$freq <- round(mean$freq, 2)

cis <- rbind(cis, mean)
rm(mean)

alpha <- c(1, .3, .3, .3, .3, .3)
names(alpha) <- c("Mean", "2011","2016","2015", "Apr 2019", "Nov 2019")

linetype <- c("solid", "dashed", "dotted", "dotdash", "longdash", "twodash") 
names(alpha) <- c("Mean", "2011","2016","2015", "Apr 2019", "Nov 2019")

cis |> 
  ggplot(aes(x = age_re, y = freq, group = year, linetype = year) ) +
  #geom_point(aes(colour = voto), size = 2, show.legend = T) +
  geom_line(aes(colour = vote, alpha= year), linewidth = 1.25, show.legend = T) +
  labs(x = "Age groups", y = "% vote") +
  facet_wrap(~ vote) +
  scale_color_manual(values = c("#2b8cbe","#d7301f", "#fe9929","#8c6bb1", "#41ae76") ) +
  scale_alpha_manual(values = alpha ) +
  scale_linetype_manual(values = linetype) +
  scale_y_continuous(minor_breaks = NULL) +
  theme_minimal() +
  theme(strip.text = element_text(face = "bold", size = 18),
        legend.position = "bottom",
        legend.title = element_text(size = 16),
        legend.text = element_text(size = 16),
        axis.text.x = element_text(angle = 45, vjust = 0),
        axis.title.x = element_text(vjust = -5),
        axis.title.y = element_text(size = 16)) +
  guides(colour = "none", 
         alpha = "none",
         linetype = guide_legend(title="Year"))

ggsave(filename = "replication/outcomes/images/FigureB6.png", height = 6, width = 9)

# C. Alternative models -------------------------------

# Figure C1
rm(list = ls())

loadRData <- function(fileName){
  #loads an RData file, and returns it
  load(fileName)
  get(ls()[ls() != "fileName"])
}

df <- loadRData("replication/database/short_time.RData")
shp <- st_read("replication/database/shp/mapa_municipios_ESPCAN.shp")
colnames(shp)[10] <- "mun_code"

map <- left_join(df, shp, by = "mun_code")
rm(df, shp)

map <- map |>  
  filter(year == "2011" | year == "2015" | year == "2016" | year == "2019")

map$depo_cat_t10[map$depo_cat_t10 == "1_no_change"] <- "No change"
map$depo_cat_t10[map$depo_cat_t10 == "2_decrease"] <- "Decrease"
map$depo_cat_t10[map$depo_cat_t10 == "3_increase"] <- "Increase"

colour <-  c("#909497","#000000","#E5E7E9")
names(colour) <- c("No change", "Decrease", "Increase")

map |> 
  ggplot() + 
  geom_sf(aes(fill = depo_cat_t10, geometry = geometry), size = NA, colour = "transparent") + 
  coord_sf(datum = NA) +
  facet_wrap(~year, ncol = 2) +
  scale_fill_manual(values = colour, name = "Population change") +
  theme_minimal() +
  theme(strip.text = element_text(face = "bold", size = 20),
        legend.position = "bottom",
        legend.title = element_text(size = 18),
        legend.text = element_text(size = 16))

ggsave(file = "replication/outcomes/images/FigureC1.png", width = 15, height = 10)

# Table C1
rm(list = ls())

loadRData <- function(fileName){
  #loads an RData file, and returns it
  load(fileName)
  get(ls()[ls() != "fileName"])
}

df <- loadRData("replication/database/short_time.RData")

setFixest_dict(c(psoe = "PSOE", pp = "PP", podemos = "UP", cs = "Cs", 
                 nswp= "PANES",  es ="ES", vox = "VOX", 
                 "depo_cat_t102_decrease"="Decrease",
                 "depo_cat_t103_increase"="Increase",
                 agriculture_workers = '% workers agriculture', 
                 services_workers = '% workers services',
                 population_log = "(Log) Population", 
                 mun_code = "Municipality", year = "Year",
                 var_older_t10 = "$\\triangle{}$% 60 y/o or older",
                 var_younger_t10 = "$\\triangle{}$% 16 y/o or younger",
                 unemployed = "% unemployed"))

m1 <- feols(c(psoe, pp, podemos, cs, nswp, es, vox ) ~ 
              sw(depo_cat_t10) +
              agriculture_workers + services_workers + unemployed +  var_older_t10 + 
              var_younger_t10 + population_log | mun_code + year,  
            data = df,
            cluster = c("region", "mun_code" , "year")) 

etable(m1, style.df = style.df(depvar.title = "", fixef.title = "", 
                               fixef.suffix = " FE", yesNo = "Yes"), tex = TRUE, signif.code = c("***"=0.01, "**"=0.05, "*"=0.10),
       order = c("Depopulation"), view = T,  replace = T,file = "replication/outcomes/tables/TableC1.tex")


# Figure C2
rm(list = ls())


loadRData <- function(fileName){
  #loads an RData file, and returns it
  load(fileName)
  get(ls()[ls() != "fileName"])
}

df <- loadRData("replication/database/short_time.RData")

modelshealth <- list()
for (i in c("psoe", "pp", "podemos", "cs", "nswp", "es", "vox")){
  modelshealth[[i]] <- df |> 
  filter(depo_cat_te %in% c("1_no_change", "2_decrease") ) |> 
    feols(as.formula(paste0(i, " ~ depo_cat_te*health +     
                         agriculture_workers + services_workers + unemployed + var_older_te + 
                         var_younger_te + population_log  |  mun_code + year")),
          cluster = c( "mun_code" ))
}

modelshealth <- lapply(modelshealth, function(x) plot_cme(x, variables = "depo_cat_te", 
                                                          condition = c( "health")) + 
                         ylab("Marginal effects") + xlab("") +
                         ggtitle(x$call$formula[[2]]) + geom_hline(yintercept = 0, linetype="dashed", color = "blue") + theme_bw())

psoe <-  modelshealth$psoe + ggtitle("PSOE")
pp <- modelshealth$pp + ggtitle("PP")
podemos <- modelshealth$podemos + ggtitle("Podemos")
cs <- modelshealth$cs + ggtitle("Cs") 
nswp <- modelshealth$nswp + ggtitle("NSWP")
es <- modelshealth$es + ggtitle("ES")
vox <-  modelshealth$vox + ggtitle("Vox")

bottom <- textGrob("\u0394 health centres", rot = 0, gp = gpar(fontsize = 20))
g_health <- grid.arrange(psoe, pp, podemos, cs, nswp, es, vox, ncol = 3, nrow = 3, bottom = bottom)

ggsave(g_health, file = "replication/outcomes/images/FigureC2.png", width = 15, height = 10)

# Figure C3

modelsedu <- list()
for (i in c("psoe", "pp", "podemos", "cs", "nswp", "es", "vox")){
  modelsedu[[i]] <- df |> 
    filter(depo_cat_te %in% c("1_no_change", "2_decrease") ) |>
    feols(as.formula(paste0(i, " ~ depo_cat_te*education +     
                         agriculture_workers + services_workers + unemployed + var_older_te + 
                         var_younger_te + population_log  |  mun_code + year")),
          cluster = c( "mun_code" ))
}

modelsedu <- lapply(modelsedu, function(x) plot_cme(x, variables = "depo_cat_te", 
                                                    condition = c( "education")) + 
                      ylab("Marginal effects") + xlab("") +
                      ggtitle(x$call$formula[[2]]) + geom_hline(yintercept = 0, linetype="dashed", color = "blue") + theme_bw())

psoe <-  modelsedu$psoe + ggtitle("PSOE")
pp <- modelsedu$pp + ggtitle("PP")
podemos <- modelsedu$podemos + ggtitle("Podemos")
cs <- modelsedu$cs + ggtitle("Cs") 
nswp <- modelsedu$nswp + ggtitle("NSWP")
es <- modelsedu$es + ggtitle("ES")
vox <-  modelsedu$vox + ggtitle("Vox")

bottom <- textGrob("\u0394 education centres", rot = 0, gp = gpar(fontsize = 20))
g_education <- grid.arrange(psoe, pp, podemos, cs, nswp, es, vox, ncol = 3, nrow = 3, bottom = bottom)

ggsave(g_education, file = "replication/outcomes/images/FigureC3.png", width = 15, height = 10)

# Table C2
rm(list = ls())
library("fixest")

loadRData <- function(fileName){
  #loads an RData file, and returns it
  load(fileName)
  get(ls()[ls() != "fileName"])
}

df <- loadRData("replication/database/short_time.RData")

setFixest_dict(c(turnout = "Turnout",
                 "depo_cat_te2_decrease"="Decrease",
                 "depo_cat_te3_increase"="Increase",
                 agriculture_workers = '% workers agriculture', 
                 services_workers = '% workers services',
                 population_log = "(Log) Population", 
                 mun_code = "Municipality", year = "Year",
                 var_older_te = "$\\triangle{}$% 60 y/o or older",
                 var_younger_te = "$\\triangle{}$% 16 y/o or younger",
                 unemployed = "% unemployed"))

m1 <- feols(c(turnout) ~ 
              sw(depo_cat_te) +
              agriculture_workers + services_workers + unemployed + 
              var_older_te + var_younger_te + population_log | mun_code + year,  
            data = df,
            cluster = c("region", "mun_code" , "year")) 

etable(m1,style.df = style.df(depvar.title = "", fixef.title = "", 
                              fixef.suffix = " FE", yesNo = "Yes"), tex = TRUE, signif.code = c("***"=0.01, "**"=0.05, "*"=0.10),
       order = c("Depopulation"), view = T,  replace = T,file = "replication/outcomes/tables/TableC2.tex")

# Figure C4
rm(list = ls())

loadRData <- function(fileName){
  #loads an RData file, and returns it
  load(fileName)
  get(ls()[ls() != "fileName"])
}

df <- loadRData("replication/database/short_time.RData")

modelslog <- list()
for (i in c("turnout")){
  modelslog[[i]] <- df |> 
    filter(depo_cat_te %in% c("1_no_change", "2_decrease") ) |> 
    feols(as.formula(paste0(i, " ~depo_cat_te*population_log +     
                         agriculture_workers + services_workers + unemployed + var_older_te + 
                         var_younger_te | mun_code + year")),
          cluster = c("mun_code"))
}

modelslog <- lapply(modelslog, function(x) plot_cme(x, variables = "depo_cat_te", 
                                                    condition = c("population_log")) + 
                      geom_rug(data = df, aes(x = population_log),  alpha = 0.2) +
                      ylab("Marginal effects") + xlab(" ") +
                      ggtitle(x$call$formula[[2]]) + geom_hline(yintercept = 0, 
                                                                linetype="dashed", color = "blue") + theme_bw())

turnout <-  modelslog$turnout + ggtitle("Turnout")

bottom <- textGrob("(Log) Population", rot = 0, gp = gpar(fontsize = 20))
g_sizemun <- grid.arrange(turnout, bottom = bottom)

ggsave(g_sizemun, file = "replication/outcomes/images/FigureC4.png", width = 15, height = 10)

# Figure C5

modelsage <- list()
for (i in c("turnout")){
  modelsage[[i]] <- df |> 
    filter(depo_cat_te %in% c("1_no_change", "2_decrease") ) |> 
    feols(as.formula(paste0(i, " ~ depo_cat_te*var_older_te + 
                        agriculture_workers + services_workers + unemployed + population_log + 
                         var_younger_te | mun_code + year")),
          cluster = c( "mun_code" ))
}

modelsage <- lapply(modelsage, function(x) plot_cme(x, variables = "depo_cat_te", 
                                                    condition = c("var_older_te")) + 
                      geom_rug(data = df, aes(x = var_older_te),  alpha = 0.2) +
                      ylab("Marginal effects") + xlab("") +
                      ggtitle(x$call$formula[[2]]) + geom_hline(yintercept = 0, 
                                                                linetype="dashed", color = "blue") + theme_bw())

turnout <-  modelsage$turnout + ggtitle("Turnout")

bottom <- textGrob("\u0394 % 60 y/o or older", rot = 0, gp = gpar(fontsize = 20))
g_old <- grid.arrange(turnout, bottom = bottom)

ggsave(g_old, file = "replication/outcomes/images/FigureC5.png", width = 15, height = 10)

# Figure C6

modelspublic <- list()
for (i in c("turnout")){
  modelspublic[[i]] <- df |> 
    filter(depo_cat_te %in% c("1_no_change", "2_decrease") ) |> 
    feols(as.formula(paste0(i, " ~ depo_cat_te*public_services +     
                         agriculture_workers + services_workers + unemployed + var_older_te + 
                         var_younger_te + population_log  |  mun_code + year")),
          cluster = c( "mun_code" ))
}

 modelspublic <- lapply(modelspublic, function(x) plot_cme(x, variables = "depo_cat_te", 
                                                          condition = c( "public_services")) + 
                         ylab("Marginal effects") + xlab("") +
                         ggtitle(x$call$formula[[2]]) + geom_hline(yintercept = 0, linetype="dashed", color = "blue") + theme_bw())

turnout <-  modelspublic$turnout + ggtitle("Turnout")

bottom <- textGrob("\u0394 Public services (education and health centres)", rot = 0, gp = gpar(fontsize = 20))
g_publicservices <- grid.arrange(turnout, bottom = bottom)

ggsave(g_publicservices, file = "replication/outcomes/images/FigureC6.png", width = 15, height = 10)

# Figure C7

models_services <- list()
for (i in c("turnout")){
  models_services[[i]] <-  df |> 
    filter(depo_cat_te %in% c("1_no_change", "2_decrease") ) |> 
    feols(as.formula(paste0(i, " ~ depo_cat_te*amenities +     
                         agriculture_workers + services_workers + unemployed + var_older_te + 
                         var_younger_te + population_log  |  mun_code + year")),
          cluster = c( "mun_code" ))
}

models_services <- lapply(models_services, function(x) plot_cme(x, variables = "depo_cat_te",
                                                                condition = c("amenities")) + ylab("Marginal effects") 
                          + xlab("") +
                            ggtitle(x$call$formula[[2]]) + geom_hline(yintercept = 0, linetype="dashed", color = "blue") + theme_bw())

turnout <-  models_services$turnout + ggtitle("Turnout")

bottom <- textGrob("\u0394 in the presence of amenities indicator", rot = 0, gp = gpar(fontsize = 20))
g_urbanfabric <- grid.arrange(turnout, bottom = bottom)

ggsave(g_urbanfabric, file = "replication/outcomes/images/FigureC7.png", width = 15, height = 10)


# D. Depopulation across region and over time -------------------------------

# Figure D1
rm(list = ls())


loadRData <- function(fileName){
  #loads an RData file, and returns it
  load(fileName)
  get(ls()[ls() != "fileName"])
}

df <- loadRData("replication/database/short_time.RData")

modelsccaa <- list()
for (i in c("psoe", "pp", "podemos", "cs", "nswp", "es", "vox")){
  modelsccaa[[i]] <- df |> 
    filter(depo_cat_te %in% c("1_no_change", "2_decrease") ) |> 
    feols(as.formula(paste0(i, " ~depo_cat_te*region +     
                         agriculture_workers + services_workers + unemployed + var_older_te + 
                         var_younger_te + population_log | mun_code + year")),
          cluster = c("mun_code"))
}

modelsccaa <- lapply(modelsccaa, function(x) plot_cme(x, variables = "depo_cat_te", 
                                                      condition = c("region")) + ylab("Marginal effects") + 
                       xlab("") +
                       coord_flip() +
                       ggtitle(x$call$formula[[2]]) + geom_hline(yintercept = 0, 
                                                                 linetype="dashed", color = "blue") + theme_bw())

g_ccaa_short <- modelsccaa$psoe + ggtitle("PSOE") + modelsccaa$pp + ggtitle("PP") + modelsccaa$podemos + 
  ggtitle("Podemos") + modelsccaa$cs + ggtitle("Cs") + modelsccaa$nswp + ggtitle("NSWP") + modelsccaa$es + ggtitle("ES") + modelsccaa$vox + ggtitle("Vox")

ggsave(g_ccaa_short, file = "replication/outcomes/images/FigureD1.png", width = 15, height = 10)

# Figure D2
rm(list = ls())

loadRData <- function(fileName){
  #loads an RData file, and returns it
  load(fileName)
  get(ls()[ls() != "fileName"])
}

df <- loadRData("replication/database/long_time.RData")

modelsyear <- list()
for (i in c("psoe", "pp", "podemos", "cs", "nswp", "es", "vox")){
  modelsyear[[i]] <- df |> 
    filter(depo_cat_te %in% c("1_no_change", "2_decrease") ) |> 
    feols(as.formula(paste0(i, " ~ depo_cat_te*year")),
          cluster = c( "mun_code" ))
}

modelsyear <- lapply(modelsyear, function(x) plot_cme(x, variables = "depo_cat_te", 
                                                      condition = c("year")) + ylab(" ") + 
                       xlab("") +
                       ggtitle(x$call$formula[[2]]) + geom_hline(yintercept = 0, 
                                                                 linetype="dashed", color = "blue") + 
                       scale_x_discrete(limits = c("1989", "1993", "1996", "2000", "2004", "2008", "2011", "2015", "2016", "2019", "2020") ) +
                       theme_bw())

psoe <-  modelsyear$psoe + ggtitle("PSOE")
pp <- modelsyear$pp + ggtitle("PP")
podemos <- modelsyear$podemos + ggtitle("Podemos")
nswp <- modelsyear$nswp + ggtitle("NSWP")
cs <- modelsyear$cs + ggtitle("Cs") + 
  geom_rect(data = NULL, aes(xmin = stage("1989", after_scale = -Inf),
                             xmax = stage("2008", after_scale = xmax-0.25),
                             ymin = -Inf,
                             ymax = Inf), color="grey", fill = "grey", alpha = .4) +
  geom_rect(data = NULL, aes(xmin = stage("2008", after_scale = xmin+0.25),
                             xmax = stage("2015", after_scale = xmax-0.25),
                             ymin = -Inf,
                             ymax = Inf), color="grey", fill = "grey", alpha = .4)

es <- modelsyear$es + ggtitle("ES") +
  geom_rect(data = NULL, aes(xmin = stage("1989", after_scale = -Inf),
                             xmax = stage("2000", after_scale = xmax-0.25),
                             ymin = -Inf,
                             ymax = Inf), color="grey", fill = "grey", alpha = .4) +
  geom_rect(data = NULL, aes(xmin = stage("2000", after_scale = xmin+0.25),
                             xmax = stage("2008", after_scale = xmax-0.25),
                             ymin = -Inf,
                             ymax = Inf), color="grey", fill = "grey", alpha = .4) +
  geom_rect(data = NULL, aes(xmin = stage("2008", after_scale = xmin+0.25),
                             xmax = stage("2015", after_scale = xmax-0.25),
                             ymin = -Inf,
                             ymax = Inf), color="grey", fill = "grey", alpha = .4) +
  geom_rect(data = NULL, aes(xmin = stage("2016", after_scale = xmin+0.25),
                             xmax = stage("2020", after_scale = xmax-0.25),
                             ymin = -Inf,
                             ymax = Inf), color="grey", fill = "grey", alpha = .4)

vox <-  modelsyear$vox + ggtitle("Vox") +
  geom_rect(data = NULL, aes(xmin = stage("1989", after_scale = -Inf),
                             xmax = stage("2015", after_scale = xmax-0.25),
                             ymin = -Inf,
                             ymax = Inf), color="grey", fill = "grey", alpha = .4)

bottom <- textGrob("Year", rot = 0, gp = gpar(fontsize = 20))
left <- textGrob("Marginal effects", rot = 90, gp = gpar(fontsize = 20))
g_year <- grid.arrange(psoe, pp, podemos, cs, nswp, es, vox, ncol = 1, nrow = 7, bottom = bottom, left = left)

ggsave(g_year, file = "replication/outcomes/images/FigureD2.png", width = 10, height = 20)

# E. Depopulation risk -------------------------------

# Figure E1
rm(list = ls())

loadRData <- function(fileName){
  #loads an RData file, and returns it
  load(fileName)
  get(ls()[ls() != "fileName"])
}

df <- loadRData("replication/database/short_time.RData")
df <- df |> 
  filter(year %in% c(2011, 2015, 2016, 2019))

data <- df %>% 
  group_by(year, depo_risk) %>% 
  summarise(mun_risk = n()) %>% 
  mutate(perc_municip = (mun_risk/8135)*100) %>% 
  na.omit()

data$perc_municip <- round(data$perc_municip, 2)

ggplot(subset(data, depo_risk == "Yes"),
       aes(x = as.factor(year), y = perc_municip)) +
  geom_col(color = "grey", fill = "grey") +
  geom_text(aes(x = year, y = (perc_municip + 2), label = perc_municip), size = 7, fontface = "bold") +
  scale_y_continuous(minor_breaks = NULL) +
  theme_minimal() +
  labs(x = "", y = "% municipalities") +
  theme(axis.text.y = element_text(size = 14),
        axis.text.x = element_text(angle = 0, vjust = 0, size = 16),
        axis.title.y = element_text(size = 16))

ggsave(filename = "replication/outcomes/images/FigureE1.png", height = 6, width = 9)
